  SUBROUTINE TRANSPOSE_XZ_XY(F,nxpsi,nxpse,nypsi,nypse,nzpsi,nzpse, &
&                     Q,imin,imax,jmin,jmax,kmin,kmax,nx,ny,nz, &
&                     comm3d,nproc,ndim,dims,coords, &
&                     nxpi,nxpe,nypi,nype,nzpi,nzpe,jsize)

! This subroutine transposes the XZ array decomposition into
! YZ decomposition

!  use MPI
  use MPI_VARIABLES, only : &
  & transpose_xz_xy_indic
  use TIMING_MOD, only : &
  & proc_time

  implicit none

  include 'mpif.h'


! Input variables       
  real(8), intent(in) :: F(nxpsi:nxpse,nypsi:nypse,nzpsi:nzpse)
  integer, intent(in) :: nxpi(0:nproc-1) ! The boundaries of the portion of array being exchanged                                    
  integer, intent(in) :: nxpe(0:nproc-1)
  integer, intent(in) :: nypi(0:nproc-1)
  integer, intent(in) :: nype(0:nproc-1)
  integer, intent(in) :: nzpi(0:nproc-1)
  integer, intent(in) :: nzpe(0:nproc-1)
  integer, intent(in) :: nxpsi,nxpse,nypsi,nypse,nzpsi,nzpse
  integer, intent(in) :: imin, imax
  integer, intent(in) :: jmin, jmax
  integer, intent(in) :: kmin, kmax 
  integer, intent(in) :: ndim, nproc
  integer, intent(in) :: comm3d
  integer, intent(in) :: nx, ny, nz
  integer, intent(in) :: jsize(0:nproc-1) ! Corresponds to basic y-decomposition of the exchanged array           
  integer, intent(in) :: coords(1:3)
  integer, intent(in) :: dims(1:3)

! Output variables
  real(8), intent(out) :: Q(imin:imax,jmin:jmax,kmin:kmax)

!      Local variables  
  integer :: i, j, k
  integer :: ki, icoord
  integer :: coordn(1:3)
  integer :: coordrecv(1:3)
  integer :: ibuf,kbuf,tag,rrank,ierr,max_maxbuf,myrank,rankrecv
  integer :: status(MPI_STATUS_SIZE,2*nproc)
  integer :: status1(MPI_STATUS_SIZE)

  integer, allocatable, save :: req2(:)
  integer, allocatable, save :: maxbuf(:) ! The maximal size of send buffer for each process 

  real(8), allocatable, save :: sarr(:,:) ! Send array
  real(8), allocatable, save :: rarr(:,:) ! Receive array
  real(8) :: pr_time, pro_time

  logical, save :: firstcall = .true.

!  pr_time = MPI_WTIME()
  call cpu_time(pr_time)

  call MPI_COMM_RANK(comm3d, myrank, ierr)

  if (firstcall) then
    allocate(maxbuf(0:nproc-1))
    allocate(req2(1:2*nproc))
    call MAXBUFCALC()
    max_maxbuf = 1
    do rrank = 0, nproc - 1
      max_maxbuf = max(max_maxbuf,maxbuf(rrank))
    enddo
    allocate (sarr(1:maxbuf(myrank), 0:nproc-1) )
    allocate (rarr(1:max_maxbuf,     0:nproc-1) )   
  endif


  tag = 1
  coordn(1) = coords(1) ; coordrecv(1) = coords(1)
  coordn(3) = coords(3) ; coordrecv(3) = coords(3)

  ki = nypi(myrank)
  if (coords(2) /= 0) then
    do i = 0, coords(2) - 1
      ki = ki + jsize(i)
    enddo
  endif

  coordrecv(2) = coords(2)
  do icoord = coords(2), dims(2) - 1
    coordn(2) = icoord
    call MPI_CART_RANK(COMM3D,coordn,rrank,ierr)
    call MPI_CART_RANK(COMM3D,coordrecv,rankrecv,ierr)
    ibuf = 1
    do k = nzpi(myrank), nzpe(myrank)
    do j = ki, ki + jsize(icoord) - 1 
    do i = nxpi(myrank), nxpe(myrank) ! Sending for coordn(2) = icoord process the jsize(icoord) portion
      sarr(ibuf,rrank) = F(i,j,k)
      ibuf = ibuf + 1
    enddo
    enddo
    enddo


    call MPI_ISEND(sarr(1,rrank),maxbuf(myrank),MPI_DOUBLE_PRECISION, &
    &              rrank,tag,comm3d,req2(icoord+nproc+1),ierr)

    call MPI_IRECV(rarr(1,rankrecv),maxbuf(rankrecv),MPI_DOUBLE_PRECISION, &
    &              rankrecv,tag,comm3d,req2(coordrecv(2)+1),ierr)


    ki = ki + jsize(icoord)
    coordrecv(2) = coordrecv(2) - 1
    if (coordrecv(2) == -1) coordrecv(2) = dims(2) - 1
  enddo


  if (coords(2) /= 0) then

  ki = nypi(myrank)
  do icoord = 0, coords(2) - 1
    coordn(2) = icoord
    call MPI_CART_RANK(COMM3D,coordn,rrank,ierr)
    call MPI_CART_RANK(COMM3D,coordrecv,rankrecv,ierr)
    ibuf = 1
    do k = nzpi(myrank), nzpe(myrank)
    do j = ki, ki + jsize(icoord) - 1 
    do i = nxpi(myrank), nxpe(myrank) ! Sending for coordn(2) = icoord process the jsize(icoord) portion
      sarr(ibuf,rrank) = F(i,j,k)
      ibuf = ibuf + 1
    enddo
    enddo
    enddo


    call MPI_ISEND(sarr(1,rrank),maxbuf(myrank),MPI_DOUBLE_PRECISION, &
    &                rrank,tag,comm3d,req2(icoord+nproc+1),ierr)

    call MPI_IRECV(rarr(1,rankrecv),maxbuf(rankrecv),MPI_DOUBLE_PRECISION, &
    &               rankrecv,tag,comm3d,req2(coordrecv(2)+1),ierr)


    ki = ki + jsize(icoord)
    coordrecv(2) = coordrecv(2) - 1
    if (coordrecv(2) == -1) coordrecv(2) = dims(2) - 1
  enddo

  endif

  call MPI_WAITALL(dims(2), req2(1), status(1,1), ierr)


  coordn(1) = coords(1)
  do icoord = 0, dims(2) - 1
    coordn(2) = icoord
    call MPI_CART_RANK(COMM3D,coordn,rrank,ierr)
    ibuf = 1
    do k = nzpi(rrank), nzpe(rrank)
      do j = jmin, jmin + jsize(coords(2)) - 1
        do i = nxpi(rrank), nxpe(rrank)
          Q(i,j,k) = rarr(ibuf,rrank)
          ibuf = ibuf + 1
        enddo
      enddo
    enddo
  enddo

  call MPI_WAITALL(dims(2), req2(nproc+1), status(1,1), ierr)

  call MPI_BARRIER(comm3d, ierr)

!  proc_time(transpose_xz_xy_indic,2) = proc_time(transpose_xz_xy_indic,2) + &
!  & MPI_WTIME() - pr_time

  call cpu_time(pro_time)
  proc_time(transpose_xz_xy_indic,2) = proc_time(transpose_xz_xy_indic,2) + &
  & pro_time - pr_time

  if (firstcall) firstcall = .false.
  RETURN
  contains

 SUBROUTINE MAXBUFCALC

 implicit none

 k = 1
 do icoord = 0, dims(2) - 1
   k = max(k, jsize(icoord))
 enddo

 do j = 0, dims(2) - 1
   do i = 0, dims(1) - 1
     coordn(1) = i
     coordn(2) = j
     coordn(3) = 0
     call MPI_CART_RANK(COMM3D,coordn,rrank,ierr)
     maxbuf(rrank) = (nzpe(rrank) - nzpi(rrank) + 1) * &
     &               (nxpe(rrank) - nxpi(rrank) + 1) * k
   enddo
 enddo

 END SUBROUTINE MAXBUFCALC
 END SUBROUTINE TRANSPOSE_XZ_XY
